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Abstract. Cultivating oleaginous microalgae in specific culturing devices such as raceways is 
seen as a future way to produce biofuel. The complexity of this process coupling non linear 
biological activity to hydrodynamics makes the optimization problem very delicate. The large 
amount of parameters to be taken into account paves the way for a useful mathematical mod- 
eling. Due to the heterogeneity of raceways along the depth dimension regarding temperature, 
light intensity or nutrients availability, we adopt a multilayer approach for hydrodynamics and 
biology. For free surface hydrodynamics, we use a multilayer Saint- Venant model that allows 
mass exchanges, forced by a simplified representation of the paddlewheel. Then, starting from an 
improved Droop model that includes light effect on algae growth, we derive a similar multilayer 
system for the biological part. A kinetic interpretation of the whole system results in an efficient 
numerical scheme. We show through numerical simulations in two dimensions that our approach 
is capable of discriminating between situations of mixed water or calm and heterogeneous pond. 
Moreover, we exhibit that a posteriori treatment of our velocity fields can provide lagrangian 
trajectories which are of great interest to assess the actual light pattern perceived by the algal 
cells and therefore understand its impact on the photosynthesis process. 



1. Introduction 

Recently, biofuel production from microalgae has proved to have a high potential for biofuel 
production [11, 36]. Several studies have demonstrated that some microalgae species could store 
more than 50% of their dry weight in lipids under certain conditions of nitrogen deprivency [11, .30, 
37] leading to productivities in a range of order larger than terrestrial plants. In this article, we focus 
on microalgae cultivation in raceways (also called high rate ponds) , whose hydrodynamics has been 
less studied than the photobioreactor culturing devices [27, 24, 29, 28, 31]. These annular shaped 
ponds of low depth (10 to 50 cm) are mixed with a paddlewheel (see Fig. 3). Due to their inherent 
nonlinear and instationary properties, where hydrodynamics and biology are strongly coupled, 
managing and optimizing such processes is very tricky. Carrying out experiments on raceways is 
both expensive and time consuming. A model is thus a key tool to help in the optimal design of the 
process but also in its operation. The objective of this paper is to propose a new model describing 
the coupling between hydrodynamics and biology within a raceway. 

The first dynamic model of a microalgal raceway pond was proposed by [-U] assuming spatial 
homogeneity. The model was later consolidated by including time-discrete photoacclimation dy- 
namics [ ]. In parallel, other less elaborated models where proposed [10, 17]. Latter on, the 
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coupling of biology with hydrodynamics in raceways was studied [ i j. ], in order to optimize the 
raceway design. In [21], algae growth and transport is modelled and several tests are performed 
in order to study for instance the effect of the water height or temperature control on the algae 
concentration evolution. However, those studies might not guarantee some key properties such as 
mass balance. We claim that the model we develop satisfies crucial mathematical properties such 
as the conservation of biological variables or positivity of the water height. 

Our approach is the following. For the hydrodynamical part, a multilayer Saint- Venant system 
with variable density introduced in [4] is applied. We use this model because we do not want to 
tackle the whole Navier-Stokes system for such shallow water situations, but we want to represent 
heterogeneity throughout the water column. This model has now been widely validated [4, 5, 32]. 
We add for our problem a specific forcing mimicking the effect of a paddlewheel. For the growth 
of microalgae, we utilize an improved version of the Droop model [12, 13, (>]. The Droop model 
has been widely studied and validated [33, 23, 7, 8]. It states that growth do not depend on the 
external nutrient concentration but on the internal cell quota of nutrients. The algae can indeed 
still grow a few days after exhaustion of the substrate thanks to their capacity to store nutrients. 
The enhanced Droop model [ti] also takes into account the effect of light on phytoplankton. We 
then write the multi-layer version of the biological model, inspired from [4]. Afterwards, we give a 
kinetic interpretation of the whole system, allowing the derivation of a numerical scheme that has 
requested properties. 

The outline is organized as follows. In section 2 we describe and justify our system of partial 
differential equations. Afterwards, we derive the numerical scheme that we will use, based on a 
kinetic interpretation. We basically follow [4] and add the new state variables concerning biology. 
In section 4 we explain how we model the raceway with the 2D code (the two dimensions being 
(x,z), respectively the length of the pool and the water depth) by adding periodic conditions. We 
also focus on the way agitation is introduced: we add a force mimicking a paddle-wheel and then we 
derive its contribution in the multilayer system, in the kinetic scheme and in the discrete scheme. 
We show that we have relevant results. Eventually, a last section deals with a Lagrangian approach 
of the algae tracking that could be useful regarding the elaboration of a better environment mod- 
eling for the algae. 



2. The coupled model 

We adopt and couple two continuous models, one for the representation of free surface hydrody- 
namics and the other for microalgae growth. We point out that it is a one way coupling: the biology 
is indeed advected and diffused by the water flow, but there is no retroaction of algae on the fluid. 
This is justified by the fact that the biological concentrations, even though greater in a raceway 
than they could be in an ocean or a lake, remain still much smaller than what is expected to change 
density or temperature of the pool (let us recall for instance that the salinity of the seawater is 
around 37 g.L~^ whereas we will never reach algae concentrations more than 1 or 2 g.L~^). 

2.1. The hydrodynamics model. Let us introduce first the hydrodynamics model in two dimen- 
sions. It will represent a free surface flow set into motion by a paddlewheel. The paddlewheel 
applies a volumic force on the water: 
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More details are given in section 4 about the formulation of the force. We begin with the 2D 
hydrostatic Navier-Stokes equations with varying density, on which we plug the paddlewheel force: 

ot OX oz ox ox oz 

i--P9+^ + ^+F.ix,z^tl (2.3) 

and we consider solutions of the equations for t > Iq, a; e M, zi,{x) < z < ri{x,t), where 
represents the free surface elevation, u = (u, w)'^ the velocity vector, p{x, z, t) is the pressure, g the 
gravity acceleration and p{T) is the water density, depending on an advected and diffused tracer T 
basically representing the temperature and which satisfies the advection diffusion equation: 

dpT dpuT dpwT d^T d^T 

The flow height is H ~ ij ~ zi,. The chosen form of the viscosity tensor is 

_ du _ du _ dw _ du 

^xx — ^M"?^~7 ^xz — ^zz — ^/^T^j ^zx — 

ox oz OZ OZ 

where p is a dynamic viscosity. 

Boundary conditions. The system (2.1)-(2.3) is completed with boundary conditions. The outward 
and upward unit normals to the free surface and to the bottom nf, are given by 

1 / _<5!Z \ 1 / 

Let Et be the total stress tensor with 

St = -pid - 

Free surface conditions. At the free surface we have the kinematic boundary condition 

|+..g-^. = 0, (2.5) 

where the subscript s denotes the value of the considered quantity at the free surface. Provided 
that the air viscosity is negligible, the continuity of stresses at the free boundary implies 

Sth, = -p'^Us, (2.6) 

where = p°'{x, t) is a given function corresponding to the atmospheric pressure. In the following, 
we assume = 0. 

Bottom conditions. The kinematic boundary condition is a classical no-penetration condition 

dzjy 

Ub.rifc = 0, or -wb = 0. (2.7) 

For the bottom stresses we consider a wall law 

tfc.STnb = KUfc.tfc, (2.8) 
where t;, a unit vector satisfying • rih = and k is a friction coefficient. 
2.2. The biological model. 
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2.2.1. Phytoplankton. Microalgae have pigments to capture sunlight, which is turned into chemi- 
cal energy during the photosynthesis process. They consume carbon dioxide, and release oxygen. 
Phytoplankton growth depends on the availability of carbon dioxide, sunlight, and nutrients. Re- 
quired nutrients are of various types, but here we focus on inorganic nitrogen, such as nitrate, 
whose deprivency is known to stimulate lipid production [25]. We consider the combined influence 
of nitrate and light on phytoplankton growth. The nutrient limitation is taken into account by a 
Droop formulation [ ] of the growth rate. The light effect (photosynthesis and photoinhibition) is 
represented using a classical formulation from ['2{i\ embedded in the model proposed by [ ]. 

2.2.2. Droop model with photoadaptation. The Droop model represents the growth of an algal 
biomass using a nutrient of concentration (nitrate under the form NO3) in the medium. 
The concentration of particulate nitrogen (nitrogen contained in the algal biomass) is denoted . 
Note that {gC.m^^), {gN.m~^) and {gN.mT^) are three transportable quantities. In 
order to reproduce the coupling between nutrient uptake rate A and growth rate /u. Droop intro- 
duced the cell quota, q{t,x), defined as the amount of internal nutrients per biomass unit: q = 
Microalgae growth rate thus depends on the intra-cellular quota: 

Mg) = Ml-y), (2.9) 

where the constant p, denotes the hypothetical growth rate for infinite quota, and Qo is the minimum 
internal nutrient quota required for growth. 

The nitrate uptake rate is a function of the external nitrate [14]: 

A(C3).A^, (2.10) 

where is the half saturation constant and A the maximum uptake rate. 

Finally, we will take into account both respiration and mortality, which are represented by a 
constant loss of the biomass with a factor R. 

In line with [(>], we modify this classical model in order to capture light and space variations. 
Light along the raceway depth. The previous model only included nutrient limited growth. It can 
be improved by introducing a new data in the system: the light intensity, which will depend on 
the quantity of water and biomass above the algae. Light is indeed attenuated by the chlorophyll 
concentration. This concentration can be linked to nitrogen through 



with 



Chl^-i{r)C^ 



I* + kj, ' 

where fc/. is a constant, /* is the average light in the water column the day before(space and time 
average). 7(/*) is presumed constant over the day. 

Let us assume that light intensity hitting the water surface is of the form 

/o = /q""'' max(0, sin(27ri)), (2.11) 

Then the intensity at depth z can be described by 

/(z) =/oe-'^(^''^*'"\ (2.12) 
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where 

^/'(C^r,z) = /" {a-f{I*)C^{z)+b)dz, (2.13) 
Jo 

(2.14) 

and /™°^ is a constant representing the maximum hght intensity, a and b are also given constants. 

The growth rate can then be computed to take into account Hght intensity, using Peeters and 
Eilers formahsm [26, 18] 

fi{q,i)^fi ^ a-^), 

with /i, Ksi, Kii three given constants derived from dedicated experiments. Finally, a down 
regulation by the internal quota of the uptake rate must be included to avoid infinite substrate (iVOs) 
uptake in the dark 

C-^ + K3 Qi 

where Qi is the maximum achievable quota. 

Adding advection and diffusion, the biological system writes in the end 

^J^c^[-^ + ^yp{^^{qJ)C'-RC'), (2.15) 

A^c^ + + P^^^C\q)C^ - RC\ (2.16) 

/^C3P^ + ^ -PA(C3,,)C\, (2.17) 



dpC^ 


dpuC^ 


dpwC^ 


dt 


dx 


dz 


dpC^ 


dpuC^ 


dpwC^ 


dt 


dx 


dz 


dpC^ 


dpuC^ 


dpwC^ 


dt 


dx 


dz 



where q = 9rr, pc^ , p.Q'^ ^ p-c^ are diffusion coefficients and (u, w) are the fluid velocities along x 



c 

and z direction. 



3. Multilayer model, kinetic interpretation and numerical scheme 

3.1. Vertical space discretization: the multilayer model. The next step is the vertical dis- 
cretization of (2.1-2.3) and (2.15-2.17) in order to obtain a multilayer system. Such a derivation 
has already been performed in [4, 5]. But in our case we have a source term representing the 
paddlewheel effect in equations (2.2) and (2.3). Since the resulting vertically discretized equations 
are not straightforward from those two previous papers, we detail it here. However, for the sake of 
simplicity, we choose to omit the viscosity terms. Moreover, we will use the classical equation of 
state relating the density and the tracer identified here with the temperature 

p(T) = po (1 - a(T - To)2) , (3.1) 

with To = 4° C, a = 6.6310^^C^ and po — 10"^ kg.m^"^. In the range of temperatures that we 
plan to use, it is easy to check that density variations are small. Hence we use the Boussinesq 
assumption which states that those density variations are taken into account in the gravitational 
force only. In any other equation the density is supposed to be po- Eventually, we end up with the 
hydrodynamic system 

du dw ^ 



dT 


duT 


dwT 




dx 


dz 




duC^ 


dwC'^ 


dt 
dC^ 


dx 
duC^ 


dz 
dwC^ 


dt 
dC^ 


dx 
duC^ 


dz 
dwC^ 


dt 


dx 


dz 
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du dv? duw I dp 1 ^ , . ,„ „, 

-^ + ^ + ^ + -^ = -Fxix,z,t), (3.3 
dt dx dz po dx po 

'^ = -pg + F,{x,z,t), (3.4) 

(3.5) 

UL U'X UZ 

and the biological system 

p{q,I)C^ - RC^, (3.6) 
= A(C^ g)Ci - (3.7) 
^~X{C^q)C\ (3.8) 

Ub UX UZ 

With g = • 

The process to obtain the multilayer system is described below. It is basically a Galerkin ap- 
proximation of the variables followed by a vertical integration of the equations. The interval [zh, ry] 
is divided into N layers {La}ae{i....,N} of thickness laH{x,t) where each layer corresponds to 
the points satisfying z e La{x,t) = [20-1/2,^(1+1/2] with 

Za+l/2{x,t) ^ Zb{x,t) +J2"=lljH{x,t), 

ha{x,t) = Za+i/2{x,t) - Za-l/2{x,t) = laH{x,t), ae[0,...,N] 

with;, >o, Ef=i^j = i- 

Now let us consider the space P^)^ of piecewise constant functions defined by 

Kh = {l.eL„(..t)W, ae{l,...,N}}, 

where ^z£La,ix,t){z) is the characteristic function of the interval La(x,t). Using this formalism, the 
projection of u, w and T onto Vq'^ is a piecewise constant function defined by 

N 

X^{x, z, t)^Yl I[-o-i/2,..+v2] iz)Xc{x, t), (3.10) 
for X e {u, w, T). The density p = p{T) inherits a discretization from the previous relation with 

N 

p^(x,z,{z„},t) = ^\z^^^,^.z^+,/.,\{z)p{To,{x,t)). (3.11) 

We have the following result. 
Proposition 3.1. The weak formulation of Eqs. (3.2)-(3.5) onV^'j^ leads to a system of the form 

E dlgH \ - dlgHUg _ I'q io\ 

a—1 a— 1 
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I o Pa+l/2 7. Pq-1/2 

Po \ OX ' ax ' 

Fx{x,z,t)dz, (3.13) 



with 



Po 

dt ^ 9a; ^^"''^'^'^'^^ ~ Ta+i/2Ga + l/2 ^ Ta-i/2Ga-\/2, (3-14) 

aG [l,...,iV]. 

a+l/2 ,„ ^ 

Ga+1/2 — 1" ""0+1/2 ^^^0+1/2 (3.15) 

Gi/2 = Gm+1/2 = 0. (3.16) 

The definitions of Pa, Pa+i/2, *^q+i/2) 7q+i/2 ol^^ given in the following proof. 

Proof. Since the demonstration for part of this proposal can be found in [4] and [5] , we wih not 
detail it here. However, we wiU focus on the integration of the agitation terms in this multilayer 
system. 

The horizontal term. The horizontal term, in the x-projection of the momentum equation will be 
handled simply: since we plan to use an expression of the force that can be integrated analytically 
(see 4.2.2), we just add the integrated force in the right-hand-side of (3.13) 

' F,{x,z,t)dz, (3.17) 

-1/2 

ae[l,...,N]. 

The vertical term. The handling of this term will require more steps. First of all, let us recall that 
in the situation where no additional source terms are present, (3.4) is used to compute the piecewise 
continuous state variable p, which is then introduced in (3.3) (see [4]). We basically proceed the 
same way. 

From (3.4), we get: 

rv rn 
p(x,z,t)^g p{x,z\t)dz'- F,Xx,z',t)dz' (3.18) 

J z J z 

If we project it on Vq'^, we get for a z in layer a 

p{x,z,t)^gi ^ Pjhj + Paiza+i/2 - z)\ ~ f F^{x, z' ,t)dz' . (3.19) 
For the vertical integration of (3.3) we need to compute: 

. , d-x'^' = ^x + (^-20) 

a — 1/2 

Since 

Pa^irj p{x,Z,t)dz, Pa+l/2 ^ P{x,Za+l/2,t), (3.21) 



8 



O. BERNARD, A.-C. BOULANGER, M.-O. BRISTEAU, AND J. SAINTE-MARIE 



we see that the vertical agitation force has an effect inp through (3.3). Let us precise the expressions 

oipa{x,t), Pa+l/2{x,t) 'And Pa-l/2{x , t) . 

Pa{x,t) ^ — p{x,Z,t)dz 



Pah. 



1/2 

N 





/■Zc,+i/2 r 


) ha , 









Pi^n-ir / F,{x,z\t)dz'dz, (3.22) 



2 

Pa+i/2{x,t) ^ g 2Z Pj^j ^ F^{x,z' ,t)dz' , (3.23) 

j=a+l "'2^ + 1/2 

Po.-i/2{x.t)^gY,Pjhj- F,{x,z',t)dz'. (3.24) 

j=a Jza-1/2 

The velocities Ma+i/2j a = 1, iV — 1 are obtained using an upwinding with respect to the direction 
of the mass exchange: 

J Ua if Ga+1/2 > , - 

We proceed in the same way for the tracer: 

(To. ifG„+i/2>0 



[ £+1 if gI+v2 < 0. ^^-^^^ 

□ 

In Prop. 3.1 the vertical velocity w no more appears, but we can derive relations for the discrete 
layer values of this variable by performing the Galerkin approximation of the continuity equation 
(3.2) multiplied by z. This leads to 



di \ 2 ] (h:\ 2 



- haWa + 2q+i/2Gq+i/2 — ^a- 1/2Gq,_i/2 , 

(3.27) 

where the Wa-, a = 1, . . . , A^, are the components of the Galerkin approximation of w on P^)^, see 
(3.10). Since aU the quantities except Wa appearing in Eq. (3.27) are already defined by (3.12), 
(3.13), (3.15), (3.16), relation (3.27) allows obtaining the values Wa by post-processing. Note that 
we use the relation (3.27) rather than the divergence free condition for stability purposes. We refer 
the reader to [ , ] for more details. 

Proposition 3.2. Multilayer version of the Droop model with photoacclimation. The weak formu- 
lation of Eqs. (3.6)-(3.8) on Pq ^ leads to a system of the form 

dh d 

+ hJp{qa,Io.)Cl~RCl), (3.28) 

dhaC^ d 

+ ho,{\{Cl,qo.)Ci-RCl), (3.29) 
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dh d 

dt^ 9a; i^aC^Ua) = C'Q+l/2G'a+l/2 - C'q_1/2G'q-1/2 

-h^X{Clq^)Cl (3.30) 

a e [l,...,iV], 

with Qa = (jt o-nd C^^^^^^, j — 1 . . .3 defined through the same upwinding as Wq+i/2 and 
(see 3.25, 3.26). 

Proof. As for Eqs. (3.2)-(3.5), we use a Galerkin approximation of the biofogical variables on P^)^. 
The Galerkin approximation on V^'^ allows to write 

N 
a = l 

for Y e (C^jC^jC"^). Then we perform an integration of Eqs. (3.6)-(3.8) over the layer a. Let us 
do it term by term. 

= — -(y„+,/,^^-r„_,/,^^) (3.32) 

n — ■ 1 /2 

For the transport terms it yields: 

^. + 1/2 ^dUahcYa . 5z„+i/2 9z„-l/2 s 

^"'^ - ("a+l/2>'o+l/2^^ - "a-l/2>'a-l/2^^) (3.33) 

a 1/2 

—K—dz = W„+l/2i^a+l/2 - Wa_i/2Ya_i/2 (3.34) 

For the reaction terms we get: 

' ^^"\ti{qJ)C^ - RC^)dz ^ h^[^l{q^J^)Ci- RCi), (3.35) 

(A(C3,g)Ci - RC^)dz = /ia(A(C;^,q„)Ci - RCl), (3.36) 

Zq-1/2 

-A(C3,g)Xdz = -/i„(A(C3,g„)Ci), (3.37) 



where 



g = — and (7„ = 



Using the notation (3.15) and (3.16) we recover Eqs. (3.28)-(3.30). □ 

3.2. Kinetic interpretation. The kinetic approach consists in linking the behaviour of some 
macroscopic fluid systems - Euler or Navier-Stokes equations, Saint- Venant system - with Boltzmann 
type kinetic equations. Boltzmann equation was first introduced in gas dynamics. It represents the 
evolution of a density of particles in a gas. Kinetic schemes have been widely used for the resolution 
of Euler equations [22, Id]. Given the analogy between Euler and Saint- Venant equations, recent 
work has been carried out to adapt those schemes to the shallow water systems ([!]). The first step 
is the introduction of fictitious particles, the definition of a density of particles and the equation 
governing its evolution. 
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3.3. The kinetic equations. The process to obtain the kinetic interpretation of the muhilayer 
model with varying density is very similar to the one used in [4] . For that reason we will only detail 
the kinetic interpretation of the biological system. For a given layer a, a distribution function 
Mo,{x,t,^) of fictitious particles with microscopic velocity ^ is introduced to obtain a linear kinetic 
equation equivalent to the macroscopic model. 



3.4. Kinetic model derivation. Let us introduce a real function x defined on M, compactly 
supported and which have the following properties 

Xi-w) = xiw) > 
Jt^x{w) dw = J^w'^xiw) dw — 1. 

Now let us construct a density of particles Ma{x,t,^) defined by a Gibbs equilibrium: the micro- 
scopic density of particles present at time t, in the layer a, at the abscissa x and with velocity ^ 
given by 

^ MM) ^ f S.-Ua{x,t) \ 

with 

and Pa defined by (3.22). 



Cq — Pa J 



Likewise, we define ^^0+1/2(2^, t, ^) by 

Na+l/2{x,t,(,) = Ga+l/2{x,t) 6 - 1ia+l/2 (a:, t)) , (3.40) 

for Of = 0, . . . , iV and where S denotes the Dirac distribution. 

The quantities 0^+1/2, < a < N represent the mass exchanges between layers a and a + 1, 
they are defined in (3.15) and satisfy the conditions (3.16), so N1/2 and also satisfy 

Ni/2{x,t,0=NN+i/2{x,t,0=0. (3.41) 

For the Droop variables, we have the equilibria 

Ui{x,t,0 = Ciix,t)M^{x,t,0, a = l,...,iV, J = l,...,3 (3.42) 

Vl+i/2i^:*,0^Cl^^^/^{x,t)N^+y2{x,t,0, a = 0,...,N, j = l,...,3. (3.43) 

With the previous definitions we write a kinetic representation of the multilayer system without 
biological reaction terms and we have the following proposition: 

Proposition 3.3. The functions (C^ ) are strong solutions of the system (3.28) -(3.30) without re- 
action terms if and only if the set of equilibria {U^{x, t, ore- solutions of the kinetic equations 

dm . 8Ui 



^'^-VLu2 + VLu2-Qul^ (3-44) 

/2 



dt dx a+1/2 ^ a-1/2 ~ ^^c/^' 
fora^l,...,N,j^l...3 with {N^+i/2{x, t, O^V^+ini^, t, 0}a=o satisfying (3.40)-(3.43) 



The quantities Q^j = Q/jj {x,t,^) are "collision terms" equal to zero at the macroscopic level 
i.e. which satisfy for a.e. values of {x,t) 

I <3c/jc'C = 0, / CQ[/.rfe = 0, and [ ^^Q^,^d^ = Q. (3.45) 
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(3.46) 



Proof. Using the definitions (3. 39), (3. 42) and the properties of tlic function x, we have 

JR JS. 

From the definition (3.40) of Na_^_i^2 we also have 

iVa+i/2(a^,i,e)d? = G„+i/2, (3.47) 



A simple integration in ^ of the equations (3.44), always using (3.45), gives the biological equations 
(3.28)-(3.30). □ 

3.5. The numerical scheme. Adding the viscous terms to the multilayer system obtained in 
Prop. 3.2 we end up with a system of the form 



dX dF{X) 



dt 



dx 



S,{X) + S,j{X) + Smo{X), 



(3.48) 



with X = (fci, ...,k],,kl...kl,kl... k%) and fci = LHCl kl = l^HCl, kl = LHCl We 
denote F{X) the flux of the conservative part, Se{X), Svj{X) and Sbio{X) the source terms, re- 
spectively the mass transfer, the viscous and friction effects and the biological reaction terms. 

We introduce a 37V x 3A^ matrix /C(^) defined by ICij — 6ij for i,j = 1,. .. ,N with Sij the 
Kroiiecker symbol. Then, using Prop. 3.3, we can write 



X = / /c(0 u\o d^, F{x) - / e/c(0 

s,{x) = / m) v\£,) dt 



(3.49) 
(3.50) 



with W{t) = [UliO, . . . , Ul^iOV, and 



Vj e 1..3. 



We refer to [5] for the computation of Svj{X). 

To approximate the solution of (3.48) we use a finite volume framework. We assume that 
the computational domain is discretized by / nodes Xi. We denote Ci the cell of length Axi = 

— with Xi^\i2 — {xi + a;i+i)/2. For the time discretization, we denote t" = X]/c<n 

where the time steps td^ will be precised later through a CFL condition. We denote 



r 3,n 



the approximate solution at time on the cell Ci with A;^'" — laH"C^^ 
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3.6. Time splitting. For the time discretization, we apply a time splitting to the equation (3.48) 
and we write 



At" dx 



= ^e(X",X"+')+5fa«o(X"), (3.51) 

S.,j{X\X-+^) = 0. (3.52) 

The conservative part of (3.51) is computed by an explicit kinetic scheme. The mass exchange terms 
are deduced from the kinetic interpretation. The biological reaction terms are not included in the 
kinetic interpretation but simply deduced from quadrature formulas. Eventually, the viscous and 
friction terms S^j in (3.52) do not depend on the fluid density p. Thus, their vertical discretization 
and their numerical treatment do not differ from earlier works of the authors [ ]. Due to potential 
dissipative effects, a semi-implicit scheme is adopted in the second step for reasons of stability. 

3.7. Discrete kinetic equation. Starting from a piecewise constant approximation of the initial 
data, the general form of a finite volume discretization of system (3.51) is 



^i+1/2 ^i-1/2 



- Ar5"+'/' + AeS^^^, (3.53) 



where cr" = Af^ /Axi is the ratio between space and time steps and the numerical flux ^7Vi/2 
approximation of the exact flux estimated at point 

In order to find a good expression of the numerical fluxes, we need to make an incursion in the 
microscopic scale. Therefore we denote the discrete particle density at time n, cell i in the following 
manner: 

77"" / £ — u" \ 
MUO = L-^x['—;r^], with<,= 

and following (3.22) 

( n I jjn N \ i-z„ + l/2 rV 

' + E PlJ^H-] -- / F^{x,,z\ndz'dz. 

Then the equation (3.44) is discretized for each a by applying a simple upwind scheme 
5ir'(0 - - ^< i^^l^i!^ - C^^i/2(0) + 

(cr/2?(e) - er/2^ (o) , (3.54) 

where 

^^.,.+1/2 j u^n^^ if£<0. 

The quantity g^^"i^^ is not an equilibrium but if we set 

l^H-+'C'^n+i ^ /s^r'COrfe (3.55) 
Jr 

we recover the macroscopic quantities at time t""*"^. 
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3.8. Numerical flux of the finite volume scheme. In this section, we give some details for the 
computation of the fluxes introduced in the discrete equation (3.53). If we denote 

= FiX^,Xl\,) = F+{X-) + F-(xr+i), (3.56) 

following (3.49), we define 



F-(xf) = / ^icmnod^, F+{xz) = / im)u:{i)d£, (3.57) 

More precisely the expression of F^{Xi) can be written 

f+{x,) = (f+,{x,),...,f+ax,), 

F+ (X,), F+^ (X,), ^^+3 {X^), . . . , F+^ (X,)) ^ , (3.58) 



with 



F+ {Xi) = Ci^,lc.H, / (u„,, + wCo,.i)x{w) dw. (3.59) 

°' J W> 



This kinetic method is interesting because it gives a very simple and natural way to propose a 
numerical flux through the kinetic interpretation. Indeed, choosing 

1 

the integration in (3.59) can be done analytically. 

However, this method proved to be numerically diffusive. That is why in practice, following 
[4, 2], we rather introduce the upwinding in the biological equations according to the sign of the 
total mass flux. We introduce then new biological fluxes that we will use instead of (3.59) 

with 

representing the mass flux in layer a at interface z + 1/2 (see [ ]) and the interfacial quantity being 
defined in the following manner 

1 C^+i if ^ < 0. 

3.9. The source terms. We refer to [ ] for the treatment of the mass exchanges terms S'e(X", X"^^) 
For the reaction terms no difficulties occur since the biological variables and what they depend on, 
light for instance, are also projected on the same Galerkin basis. We write 

/ / i?(x,z,t") = /i„Aa;,i?^ (3.61) 

where i?(x, z, <) represents one of the three reaction terms we have in the Droop model. 
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3.10. Properties. Although we are not gomg to provide any proof in this section we want to 
recall some important features of the model used in this paper (see [4, 5, 32] for the detailed 
proofs). First of all, under a certain CFL condition which means that the quantity of water leaving 
the cell during a time step is less than the current water volume of the cell, the water height and 
the biological concentrations remain non-negative. Second of all, [4] states that a passive tracer 
satisfies a maximum principle. Finally, we want to specify that a second order scheme in space 
and time is possible and has been performed in the simulations of section 4. The second order in 
time is achieved through a classical Heun method. We apply a second order in space by a limited 
reconstruction of the variables [3] based on the prediction of the gradients in each cells, a linear 
interpolation followed by a limitation procedure. 



4. Simulations 

4.1. Analytical validation on non trivial steady states. In this section we show that we 
can find analytical solutions of the coupled problem, provided that we use Euler equations and 
a simplified model for biology. We want to emphasize the importance of this part. Validating a 
numerical code is indeed a complex but highly required task when it comes to non trivial situations. 
It is clear though that analytical solutions of the whole coupled problem cannot be found. However, 
we propose here a simplified version of a biological model that would be embedded in free surface 
Euler equations and for which we can find non trivial steady states. 
Let (u,w) be the following vector field: 

ui.,z) = a r'^^^;'^f , (4.1) 
sm(pH) 



a/3 

w(x, z) 



sin{j3{z — Z}j))cos{PH)— h cos{l3{z — Zb))sin{j3H)—— 

ox ox 



(4.2) 



with a, j3 being given constants and H{x) being known. In ['j], the authors show that it is an 
analytical solution of the Euler system of equations for free surface flows. In order to validate 
the coupling with biological equations which is performed in this paper, we propose a simplified 
stationary coupled model. 

Let us consider the following scalar field 

r(x,z)=e-(^-(---^». 

We can check by simple calculation that T is solution of 

duT dwT 



^ dx dz 

where 



= f{x,z)T{x,z) (4.3) 



1 reaction 
advection 



Clearly, this equation is a simplified and alternate version of a biological model, with advection and 
reaction terms. Notice that solutions to (4.3) represent non trivial equilibria between advection 
and reaction terms. 

Eventually, we compare analytical and numerical results of the following situation: a given hy- 
drodynamic and biological fiow is imposed at the left boundary of a 20m long pool with topography 



A 2D MODEL FOR HYDRODYNAMICS AND BIOLOGY COUPLING. 



15 



Zb = 0.2e'-^~*^ — 0.4e'-^~^^^ , a given water height is fixed at the right boundary. To solve this an- 
alytically we follow [!)]: given the topography Zf,, between horizontal coordinates x = and x = 20, 
we recover H{x) and deduce u{x,z) and w{x,z) for a = 0.4 and /3 = 1.5 thanks to the above 
formula (4.1, 4.2). Numerically, we impose the analytical flux at the left boundary, the analytical 
height at the left boundary, and run the code for 500s, with 300 nodes and 20 layers (thus 6000 
vertices), while initial conditions were set up to zero for every variable (u, w, T). We see in Fig.l 
that we recover numerically the hydrological and tracer steady states solutions for this simplified 
coupled model. Therefore, we consider that our method and our code are likely to produce valid 
results. 



u(m/s) T(g/L) 




Figure 1. Analytical (top) and numerical (bottom) steady states solutions for hy- 
drodynamics and tracer in case a = 0.4, /3 = 1.5 and Zb = 0.2e^^~*^ — 0.4e*^^~^^-' . 
A given flow is imposed at the left boundary whereas the output boundary condi- 
tion concerns the water height. 



4.2. Numerical simulations of a raceway. 

4.2.1. Light. As explained in section 2, the light intensity at a particular point will depend both on 
the amount of water above this point (i.e. related to the depth) and on the quantity of encountered 
chlorophyll (proportional to the concentration in intracellular nitrogen). Additionally, we explain 
in section 3 that a Galerkin approximation of the variables is carried out. Hence light is also 
discretized along the depth by layers. We show in Fig. 2 the light profile for different values of 
7(/*). Recall that this value is related to the average irradiance perceived by microalgae the day 
before (photoadaptation phenomenon). The parameters used are exposed in Table 1. We see for 
any curve the the exponential decay. Moreover, this figure shows that the more exposed to light 
the microalgae were, the more receptive to it they are (no photoinhibition occurs in those ranges 
of irradiance). 
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Light intensity in the water column 





0.25, 


' = 80.64 (imol.m"^.s"^ 


— 7 = 


0.20, 


= 1 16,55 fimol.m"^.s" 


— y= 


0.15, 


'= 176.40 fimol.m"^.s" 


— Y = 


0.10, 


* = 296.10 (imol.m"^.s" 




Figure 2. Light intensity 
in the water column for 
a concentration in nitrogen 
equal to 5.0 gN.m~^. The 
figure shows the exponential 
decay of light with depth for 
different values of 7(/*). 



Parameter 


Value 


Unit 


^Q,rnax 


500 


IJLmol.m~^ .s^^ 


a 


16.2 


m^.gChl-^ 


b 


0.087 






5.0 


gN.m^^ 



Table 1. Parameters used 
for the computation of light. 



4.2.2. Paddlewheel. We aim at representing the kind of raceway shown in Fig. 3. In order to model 
it in 2D, periodic conditions are applied on a rectangular pool containing a paddlewheel in its first 
half (Fig. 4). The wheel is not modelled physically but is represented by a force that is able to mimic 
its effect and give the system an equivalent energy. Therefore the following expression is assumed: 



"x{x,Z,t) = F (^y^{x~ X^heelY + {z ~ Z^heelYuj^ COS 9 (4.5) 
\{X, Z,t) ^ F (x -~ X^heelf- + {z ~ Z^heelf^ sin 6* (4.6) 



where F is a constant, Q is the angle between the blade and the vertical direction, uj ~ 6. The force 
is normal to the blade and depends on the square of the velocity of the point on the blade: the 
further it is from the center of the wheel, the bigger is the force. Therefore, the energy provided by 
the wheel is be proportional to the cube of the velocity. To include it in the model, the process is 
similar to what is explained in section 3.5: the force is added in the Navier-Stokes equations in order 
to derive again the multilayer model. This adds a source term in the x-momentum equation and 
change the expression of the pressure given by the z-momentum equation. We remind the reader 
that non-hydrostatic terms have not been taken into account in the derivation of the model (though 
it is described in [•)2]). The obtained results in Fig. 5 let us think that this may be appropriate to 
model the paddlewheel effect. 

We performed several simulations for different paddlewheel angular velocities. We use a pool with 
length 20m and height 0,5m. We show in Fig 5 that we are able to capture realistic hydrodynamics: 
a laminar flow of reasonable horizontal speed far from the wheel and a turbulent flow close to it. 
Concerning the hydrodynamics parameters, we take into account horizontal and vertical viscosity 
(/i = 0.00lTO^.s~^)and a Navier-type bottom friction( k — O.Olm.s"^). Besides, lj = 0.85rad/s. 



A 2D MODEL FOR HYDRODYNAMICS AND BIOLOGY COUPLING. 



17 




Figure 4. Left: Pool and wheel dimensions and positions. Right: Outlook of the 
force applied to model the effect of the paddlewheel, supposed to be located in the 
first half of the pool, with maximum efficiency at the end of the blade. 
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a 



Paddlewheel located cell 



Cell far from paddtewheel 




Horizontnl velocity - 0.6 tt.04 iri/s 



Mi 



u(m/s) 

w(nVs) 



Vertical vetocitv - +/- 0.05 m/s 

iiiiiiiiiifiiiiiiiiiiiiiiiiiiniiliiiii iii'inii iiiiiiB 



3lo 35 Bie 

Time(s) 



"IS 9^10 




Hori;:onliil velQc[lv u 0.48 m/s 



-u(rr/s) 



Vertical velocity w ~ m/s 



ittB <3ii 55jj iic 

Time(s) 



"135 



Figure 5. (a): Velocities along vertical and horizontal axis in a cell located near 
from the wheel rotating at angular speed w = 0.8rad/s. The flow is very turbulent, 
(b): Velocities along vertical and horizontal axis in a cell located far from the wheel. 
An asymptotic value of OASm.s^^ is reached. 



In order to have a first idea of what is happening in the fluid, we add a tracer which is advected 
and slightly diffused. Fig 6 illustrates the effect of the wheel on the mixing. In particular, after 
several minutes, the pool seems well homogenized (here again uj — 0.85rad/s). 

4.3. Results. Eventually, we performed numerical simulations of the whole coupled model ex- 
plained in section 2 and discretized in section 3. We use the same pool with length 20m and height 
0,5m. The height is chosen such that at the bottom of the pond, the respiration rate is close to 
the growth rate. We perform several 20 days simulations for different agitations, different initial 
conditions and we compare them to reference simulations i.e. without paddlewheel. The parameters 
used for the simulation (concerning the biological system) are exposed in Table 2. The relevant 
details of every simulation are in Table. 3. The results are depicted in Fig. 7. 

The initial concentrations of particulate and dissolved nitrogen are the same for the six simu- 
lations. The initial microalgal carbon concentration varies (the initial internal nitrogen quota is 
therefore changing as well), and different agitation velocities are tested. The plots represent the 
average concentrations in the raceway (see local concentration comparison between upper and bot- 
tom layer in Fig. 8). The carbon curves (Fig.7-a) show in every situation that agitation leads to 
better productivity. This is explained by the lack of nutrients in the upper layers at some point. 
However, regarding the initial internal quota (Fig.7-b), the time when agitation actually enhances 
productivity varies. This phenomena is due to the fact that the internal nutrient pool, if not filled 
enough(low quota q), will tend to increase further by absorbing more external nitrogen. Thus lead- 
ing two main consequences: the extracellular nutrient concentration diminishes and the intracellular 
nitrogen increases. Since chlorophyll is positively correlated with the latter biological variable, the 
light can not penetrate so deep anymore. We point out that the model is not able to differentiate 
between several agitation velocities. Ideas to explain and overcome this limitation are suggested in 
section 5. 

From these series of simulations we can conclude that our model is capable of reproducing 
coherent results, in the hydrodynamical part (adequate asymptotic velocity, turbulences near the 
wheel and laminar flow far from it) as well as in the biological concentrations. Nevertheless, we can 
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t = 2000 s 



le-09 



Figure 6. Snapshots of the tracer (poUutant) concentration in a raceway set into 
motion by a paddlwheel of angular velocity w — O.Sbrad/ s. It is clear that after 
several minutes, the raceway is totally homogeneous. Therefore, the paddlewheel 
has indeed the required effect on the mixing. 



only consider those results as preliminary, since no quantitative comparison with any data has been 
provided. Future work would include adaptation of the biological model given the hydrodynamical 
results obtained in this study. Moreover, other variables should be added to the system. For instance 
temperature, which could have a non negligible effect, sedimentation etc. Finally, experimental data 
should help us calibrate more accurately the raceway parameters and extend it to three dimensions. 



5. LAGRANGIAN APPROACH 

In order to better understand the physiology of the microalgae submitted to variable light and 
nutrient stresses, and improve the biological models in these realistic conditions, it is crucial to 
study these organisms at individual scale when excited by such signals whose frequency is much 
faster than in natural environments such as lakes or oceans. However, it is very tricky to measure 
the light and nutrient signal perceived by a cell in a raceway. A Lagrangian approach derived from 
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Parameter 


Value 


Unit 




1.7 




Qo 


0.050 




Qi 


0.25 


gN.gC-^ 


Ku 


295 






70 


limol.m"^ .s~^ 


A 


0.073 


gN.gC-\day-^ 


Ks 


0.0012 


gN.m~^ 


R 


0.0081 


day~^ 




500 


IJ,mol.m~^.s~^ 




0.25 


gChl.gN-^ 


a 


16.2 


m^.gChl-^ 


b 


0.087 





Table 2. Parameters for the simulations. 





Simu 1 


Simu 2 


Simu 3 


Simu 4 


Simu 5 


Simu 6 


CoH.9.m-3) 


25 


25 


50 


50 


83 


83 


{^)^i9N.gC-^) 


0.2 


0.2 


0.1 


0.1 


0.06 


0.06 


Cl{g.m^'^) 


5 


5 


5 


5 


5 


5 


Agitation 


no 


yes 


no 


yes 


no 


yes 


C)(5.m-3) 


60 


74 


79 


103 


100 


129 


(fi)^{gN.gC-') 


0.115 


0.120 


0.110 


0.085 


0.087 


0.065 


C%g.m-^) 


3.77 


0.0 


0.0 


0.0 


0.0 


0.0 



Table 3. Initial conditions and final concentrations for the set of 6 simulations 
we performed. The differences lie on the one hand on the fact that the water is 

agitated or not and on the other hand on the initial carbon concentration, which 
changes the internal quota (nitrogen concentration is always equal to 5{{g.m~^)) 



the previous model can lead to the reconstruction of cell trajectories. Prom this information, it is 
then possible to derive the signals that microalgal cells undergo. 

In order to follow the position of a particle, we simply need to integrate the following equation. 
If M{t) is the position of particle M at time t, then we have: 

^«=..(M,„,0 (5.1) 

where v{x, t) is the eulerian velocity field at position x, time t. We do not add to those particles any 
Brownian motion that would refer to the diffusion of biological concentrations at the macroscopic 
scale. Therefore, more realistic Lagrangian trajectories should look more noisy than what we get 
but the general behaviour is well represented. 
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^1^=0.06 agitation 
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Figure 7. (a): Carbon concentration; (b): Internal quota q; (c): nitrogen con- 
centration; (d): substrate concentration(A^03). Those plots illustrate the average 
concentrations in the raceway for 6 simulations. Three were carried out without 
agitation, and the other three had the agitation term. In each situation, agitation, 
leading to homogenization leads to a better productivity. However, for certain 
initial conditions, the improvement is quite slow (after several days), since the 
biological variables do not evolve as quickly as hydrodynamics does. 



The global simulation in the previous section could not differentiate between several paddlewheel 
velocities. It is likely that the biological model is not affected, provided that the mixing is efficient. 
Therefore we try to understand what could be those differences. For several angular velocities, we 
perform a one hour simulation and build the trajectories afterwards. One hundred particles are 
equally distributed in the pool along the depth dimension. Fig. 9 shows the distribution of particles 
against the percentage of time spent under more than 50% of the incident light (we will call it 
high enlightenment). We notice that the wheel velocity has an influence on the number of particles 
which never undergo high enlightenment (38 particles over 100 for uj — 0.5 and only 13 particles 
over 100 for lu = 1.0). But globally, the particles are distributed around 20% of their time under 
high enlightenment. This is what we can expect from a good mixing since the light is exponentially 
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^ Local carbon concentration (C^ for q,^., = 0.2 gC.gN 

180r <* 




'o 2 4 6 8 10 12 14 16 18 20 
Time{days) 




'o 2 4 6 8 10 12 14 16 18 20 
Time(days) 



Figure 8. (a): Local carbon concentration when qinit — 0-'^9N.gC~^ , in the bot- 
tom and upper layer, with and without agitation, (b): Local carbon concentration 
when qinit = O-OGgN.gC^^ , in the bottom and upper layer, with and without agi- 
tation. We see clearly the homogenization due to agitation. In both cases, we also 
see that the bottom layer when not agitated does not vary so much, which means 
that we are close to the point when respiration compensates growth. 



decreasing and getting 50% of the incident light means being in the first 10cm from the surface of 
the pool (over 50cm). 

In Fig. 10 we plot different indicators for eight velocities, established during a one hour simulation. 
First of all, Fig. 10. a illustrates the average velocity at the end of the simulation. It turns out to lie 
in realistic ranges (0.2 m.s~^ to 0.7 m.s~^). Second of all Fig. 10. b represents the average proportion 
of time spent for any particle under high enlightenment. From this curve we deduce that the global 
percentage of high enlightenment may not vary much between uj — 0.5 and lo — 1. The mixing is 
indeed well carried on (as shown in previous section) and in average, the particles have the same 
history. Finally, Fig. 10. c depicts the number of times the particles switch from low enlightenment 
to high enlightenment during one hour. This last indicator gives insights about the duration of the 
high enlightenment periods. The greater it is, the shorter but more numerous were the high light 
instants (the particles switches more often). 

Finally, the trajectories of three particles are depicted in Fig. 11. a. From those trajectories we 
can extract two important informations. First of all, between two passages around the paddlewheel, 
the particle seems to stay at constant depth, thus enforcing the fact that the flow is laminar apart 
from the wheel. Second of all, we clearly see that the depth of the particle is suddenly modified by 
the wheel, giving rise to abrupt changes in the enlightenment. Fig. 11. b shows the light received by 
those three particles in the case where the average intracellular nitrogen concentration is 5 gN.m~^ 
and 7(/*) — 0.1 gChl.gN~^ (high irradiance the day before the simulation). 

Fig. 2 shows that light intensity is very low in the 20 deepest centimeters. Regarding our results 
on Fig.ll, we assume that microalgae should spend one lap at high light and then several laps at 
low light. Since the asymptotic velocity is in the range 0.3m. — 0.8m. s~^, a whole lap is done in 
the order of a minute. Therefore, algae are faced to light changes with a time scale in the range of 
ten minutes. As regards our modeling problem, this tells us that a biological model valid for these 
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Figure 9. Number of particles for each class of enlightenment for four wheel 
angular velocities (0.5, 0.75, 0.85 and 1.0 rad.s~^). One class represents a range of 
5%. Being for instance in the class 20%-25% means that the particle spent between 
20 and 25% of her time under hight enlightenment (more than 50% of the incident 
light). 



fast time scale could lead to more precise results. Additional experiments forcing microalgae with 
typical light signal deduced from Figure 2 must therefore been carried out and support a microalgae 
modelling at fast time scale [15]. 

6. Conclusion 

In this paper we derived a new model coupling hydrodynamics to biology in two dimensions. It 
provides new insights to better understand and represent this nonlinear and non stationary complex 
process. The multilayer model adequately represents the high vertical heterogeneity characterizing 
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Average velocity magnitude Time spent under more than 50% of the incident iight (%) 




§.5 0.55 0.6 0.65 0.7 0.75 0.8 0,85 0.9 0.95 1 
Omega(radian.s~^) 



Figure 10. (a): Average velocity after one hour for 8 velocities, (b): Percentage 
of time spent under more than 50% of Iq. (c): Number of times a particles switches 
from a situation where it perceived less than 50% of light intensity to a situation 
where it perceives more during one hour. 

the agitated raceway. A key point of our approach, is that we were able to provide analytical 
solutions which validated its numerical integration. The results show that water agitating through 
the paddlewheel has an effect on the growth of algae, particularly because of the lack of nutrients 
at the surface versus the lack of light around the pool's floor. One of the outcome of this work 
is the identification of realistic light signals to which microalgae are faced. Lab scale experiments 
will be performed to assess the impact of such high frequency light signals on microalgae. It is 
worth remarking that similar works have been carried out for photobioreactor [27, 2(S], leading 
to the identification of much faster time scales in the range of the second. The comparison with 
experimental data would be of great relevance to better calibrate the hydrodynamics, and later on 
improve the biological predictions of the model. It is clear however that many parameters need to 
be taken into account in order to increase the model prediction capacity, for instance temperature. 
It was here considered as a passive tracer, only advected and diffused, with no particular effect on 
the biology (temperature does not appear in the growth or respiration rate). The sunlight effect 
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Figure 11. (a):Trajectories of three particles during the simulations. The large 
curve represents the water surface at the middle of the pool. The other plot is the 
height of a given particle through time. The algae undergo sudden changes of depth 
every time it meets the wheel. (b): Perceived light from the microalgae. Particles are 
subject to even greater irradiance changes since the light is exponentially decaying. 



on water temperature will be taken into account in a next stage since it deeply affects microalgae 
growth. Moreover, some microalgae species do not swim in the water and tend to sediment. This 
property could increase the beneficial effect of the wheel compared to a situation where the wheel 
is very slow or absent. Computational time is another issue which has to be improved in order 
to use the model e.g. for process optimal design. Up to now, schemes we are developing are only 
explicit. We are therefore constrained with a restrictive CFL condition. Improvement towards 
implicit schemes is also a future concern. 
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